Disordered driven coupled cavity arrays: Non-equilibrium stochastic mean-field 

theory. 
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We study the interplay of disorder with pumping and decay in coupled qubit-cavity arrays, the 
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ity present in the clean pumped system, and that moreover, the combination of disorder in on-site 
energies and decay can lead to effective phase disorder. To explore these questions, we present a non- 
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questions. This technique is developed for rather general forms of light-matter coupling, driving, 
dissipation, and on-site disorder, making it applicable to a wide range of systems. 

PACS numbers: 42.50.Pq, 72.15.Rn, 03.75.Kk 



I. INTRODUCTION 

Quantum simulation [1] concerns how controllable 
quantum systems can be used to model particular desir- 
able Hamiltonians, in order to find the ground state, or 
other properties, of otherwise hard-to-simulate problems. 
Recently, there has been significant progress in realizing 
quantum emulators based on systems including ultracold 
atoms [2| , Rydberg atoms [3| , Trapped Ions 0, Q or su- 
perconducting qubits in microwave cavities [6|-[8j. One 
approach that has been used recently for cold atoms is 
to engineer an effective Hamiltonian in a rotating frame, 
by using a Raman pumping scheme This approach 
has been used to realize the superradiance transition in 
the Dicke model [lOj. In such cases, and more gener- 
ally for coupled matter-light systems such as supercon- 
ducting qubits in microwave cavities, it can be crucially 
important to understand the effects of losses and dissi- 
pation. For example, in the Dicke model, the presence of 
losses means the critical behavior [111— Tl3[| becomes clas- 
sical [14] due to the effective temperature introduced by 
losses. Similar issues can generally be expected to oc- 
cur in any open driven system, and this therefore may 
have consequences across the range of experimental sys- 
tems considered as potential quantum emulators, and in 
particular, for coupled light-matter systems [15j. 

Another model where quantum simulation has been 
explored 1161 is the disordered Bose-Hubbard model 
(BHM) [17j. This model consists of bosonic particles 
hopping between sites with repulsion between particles 
on the same lattice site. This model can be simulated 
with ultracold bosonic atoms, introducing disorder in a 
highly controlled manner by superimposing a fine-grained 
optical speckle potential with a periodic optical lattice 
[18l-l2l|. In the presence of weak disorder in the on-site 
energies, three possible ground states exist at zero tem- 
perature: a superfluid phase and two insulating phases. 
The two insulating phases are the incompressible Mott 
Insulator and the compressible Bose glass. In the Mott 
Insulator, the particles are localized because of strong 



local repulsions, in the Bose glass particles are localized 
because of the disorder potential. Despite the long his- 
tory of the BHM, it is only recently that several aspects 
of this model have been fully understood, such as confir- 
mation that the Bose-glass phase always intervenes be- 
tween Mott insulator and superfluid [22| , and the distinc- 
tion between the Mott insulator and Bose glass regard- 
ing whether fluctuations are self- averaging [23]. Even if 
quantum simulation of such a model with an effective 
Hamiltonian in a dissipative system can only model the 
finite temperature case, this may of itself be enough to 
answer questions such as whether the finite temperature 
insulating phase is self- averaging. However, as we will 
discuss further below, disorder and non-equilibrium ef- 
fects can conspire to significantly change the behavior 
(and universality class) of the model system. 

The Jaynes-Cummings-Hubbard model (JCHM) [24| is 
very closely related [6] to the BHM, and can be more di- 
rectly realized in coupled cavity arrays [6-8,[l5|,[25|. This 
model naturally describes superconducting qubits in mi- 
crowave cavities. The JCHM consists of photons coupled 
to two- level systems, considering photons confined in an 
array of coupled cavities, with weak hopping between dif- 
ferent cavities. The JCHM requires that only the energy 
preserving term (a<r + + aV~) in qubit-cavity coupling is 
important. When counter-rotating terms (aa~ + aV + ) 
are also important the model is known as the Rabi- 
Hubbard model, and the symmetry of the Hamiltonian 
is lowered from Z7[l) to Z<i and the phase diagram signif- 
icantly changes [26|. In the case of the JCHM, previous 
work has shown how in equilibrium, including on-site dis- 
order leads to behavior very similar to the BHM [2?}, as 
may be expected from the symmetries of the problem [(| . 

In this paper we study the non-equilibrium JCHM in 
the presence of disorder. We focus on the simplest possi- 
ble form of pumping and decay, i.e. uniform coherent 
pumping, as has previously been studied in the clean 
limit [28l-l3Q( . In this case, all symmetries are broken 
by the pumping, and no phase transitions are expected. 
Nonetheless, the behavior we observe and discuss for this 
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case clearly shows how new physics would also arise with 
other forms of pumping which need not break the sym- 
metries of the model. In particular, we see that pump- 
ing and dissipation can transform on-site energy disorder 
into phase disorder, destroying long range order in the 
supernuid phase. In addition, we explore the fate of the 
bistability seen in the clean non-equilibrium JCHM (28| . 

To explore these questions, we use a generalized 
"Stochastic" Mean-Field theory [3ll-[33|, which involves 
self-consistency equations for the probability distribution 
of local order parameters. We extend this approach to 
apply to open quantum systems. Such an approach is 
approximate, and only becomes well controlled at high 
coordination number (i.e. in high dimensions). Nonethe- 
less, it provides a simple tool to effectively explore the 
interplay of disorder and pumping, and see whether ef- 
fective Hamiltonians for open systems could in principle 
be used to simulate disordered quantum systems. 

The remainder of this paper is arranged as follows. 
In Sec. HIl we generalize stochastic mean-field theory 
(SMFT) to treat disorder in open quantum systems. The 
technique is introduced for rather general forms of pump- 
ing, decay, and on-site disorder. As an example, we apply 
the SMFT to the dissipative JCHM. In Sec. EH we first 
briefly summarize the behavior of the JCHM in the ab- 
sence of disorder and then discuss the effects of on-site 
disorder in the excitation energies of the two-level sys- 
tems. 



II. STOCHASTIC MEAN FIELD THEORY OF 
OPEN SYSTEMS 

This section briefly summarizes the SMFT approach as 
applied to the non-equilibrium problem. The equilibrium 
SMFT was introduced in the context of disordered anti- 
ferromagnets [31] and later applied to the BHM [32|, l33| 
and has more recently been applied to the JCHM |27j . 
We present the following discussion for a general cou- 
pled cavity array problem, and specialize to the JCHM 
in section ITTTI 

We consider an array of cavities with coordination 
number z and hopping J/z of photons between neigh- 
boring cavities, given by the Hamiltonian 

i (ij) 

where a\ (a$) creates (annihilates) a photon on the zth 

cavity. The on-site Hamiltonian hi = h(ai,x\ a \ei) for 
the individual cavities can be completely general at this 
point. The operators X^ act on the Hilbert space of the 
possible quantum states of the matter contained in the 
cavities. In the simplest cases, including the JCHM and 
the Rabi-Hubbard model, this will be a two-level system 
and the X operators are spin-1/2 operators. The on-site 
Hamiltonian will contain a coupling between the photons 



and the matter degrees of freedom as well as any coherent 
pumping terms. 

We further introduce on-site disorder e$ which can cou- 
ple either to the photon energy or to the matter in the 
cavity. The disorder follows a probability distribution 
p(e) and is assumed to be uncorrelated between different 
cavities. For such on-site disorder the method is as devel- 
oped in Refs. [32|, |33|. If one instead considered disorder 
in the hopping between sites, the problem is analogous 
to that originally considered in Ref. [3l|. 

Dissipation is included on the level of a master equa- 
tion for the time evolution of the density operator, 

% = -i[H,p]+Y,\ ? £ N + £ Y c[x * a)] 1 ' (2) 

i I a ) 

where C[X] = 2XpX^ - {X*X,p} denote the standard 
Lindblad operators. 

The basic idea of SMFT is to consider a self- 
consistency condition for the probability distribution 
P(i/j) of on-site coherent fields ip = (a). From this one 
may find the distribution of sums of fields from neighbor- 
ing sites: 

w) = / n d ^ - E ^ wo (3) 

where the product and sum run over the z nearest neigh- 
bors. The relation between P and Q simplifies in Fourier 

space 

Q(0 = J d<PQ(4>)e^, Q{4>) = J ^Q(0e"** 

(4) 

for the Q distribution. Using the convolution theorem 
we obtain Q(f) = P{£) z . 

Given the distribution of fields from neighboring sites, 
the self-consistency condition comes from assuming that 
this distribution of fields is uncorrelated with the site 
energies, and so one may write: 

P(V>) = Jd<f>J deQ(4>)p(e)5^ - \(<f>,e)). (5) 

Here A(e, <p) gives the expectation for ip corresponding to 
a field <j) from the neighbors, and on-site energy e. In our 
case this corresponds to finding the steady-state on-site 
density matrix from, 

^ = -i [hf, Pi ] + fa] + £ ^C[X&] (6a) 

a 

hf = h i~ J - (<Ml + 0* a i) ( 6b ) 

and determining the expectation value A(0, e^) = 
Tr(diPi). In steady state, dpi/dt = 0, the master equa- 
tion for the onsite density operator (j6j) turns into a 
set of coupled linear equations for the matrix elements 
(pi)mn — (m\Pi\n) with respect to a basis of the prod- 
uct Hilbert space of the matter and photon systems. 
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While the former is usually finite, we truncate the bosonic 
Hilbert space at a certain maximum number of photons 
per cavity. 

As noted above, even for real (/>, the values of ip will be 
complex. This means it is necessary to allow for the dis- 
tributions P and Q to extend over complex fields. Convo- 
lution of two-dimensional distributions follows as before, 
but with <j) — > ((/)',(/)") in order to use the convolution 
theorem. In practice, we find it most efficient to pre- 
calculate an interpolated approximation to A(</>, e), and 
then iteratively update P(V0, Q{(j>) until the distribution 
converges. One may also note that it is not guaranteed 
that the above iteration procedure should converge, nor 
that it should converge to an unique solution — as dis- 
cussed below, the mean field decoupling introduces the 
possibility of multist ability. However, in cases where it 
does converge, the solution found can be regarded as an 
approximate description of a possible asymptotic state of 
the system. When multiple solutions exist, further work 
is required to determine which solution is reached from 
given initial conditions, and the rate of tunnelling events 
that may switch between solutions. This is discussed 
further below in Sec. IIII A[ In the cases presented in this 
paper, only one asymptotic state was found. 



III. APPLICATION TO THE PUMPED 
DISSIPATIVE JCHM 

As a simple application of the above technique, and the 
simplest kind of pumped-dissipative array, we consider 
here the coherently pumped Jaynes-Cummings-Hubbard 
model, as studied previously in [28l430| . In terms of the 
general lattice problem described in Eq. (pQ), the Jaynes- 
Cummings-Hubbard model that we consider has an on- 
site Hamiltonian: 



JCHM has a Master equation with 



hi = Ja\ai + -^-of + g[cft a i 



H.c.) 



H.c.) (7) 



where a\ creates a photon in the zth cavity and the spin- 
1/2 operators a+ , a~ describe transitions of the state of 
the two- level (artificial) atom on site i. f denotes the 
strength of the pumping at frequency uj p . The cavity 
photon energy J is chosen so that for g = 0, the bottom 
of the photon dispersion is at zero energy. Disorder is 
introduced by considering a Gaussian distribution of e^, 
of width <j e and we take the mean value e = 0, so that 
the mean detuning is as in [28]. Further, we consider loss 
terms of the form £\ {(n/2)C[ai} + (7/2)£[crr]}. The 
problem can be trivially made time-independent by the 
Unitary transform a — >• ae~ lUJpt ,a~ — >• <j~e~ lUJpt . 

Other than the coherent pumping term, the problem 
we consider has a (7(1) symmetry, and this can be used 
to simplify the pre-calculation of A(0, e) as discussed in 
the previous section. The effective on-site problem of the 



hf = ( J - u p ) a}ai + 



-a- + g(a+ai + H.c.) 



/ - Ui + H.c. 



(8) 



One may then write the steady state expectation Tr(pa^) 
arising from this effective Hamiltonian along with the 
Lindblad terms in the form: 



TrM = A I / eff = / 



J</> 



(9) 



where the last line of Eq. 8 can be written as 
^yeffj* a . _|_ J eff a jj 5 combining both the explicit pump 
and the field coming from the neighboring cavities 
into / eff . The advantage of writing the expression 
in this form is that one may note that A(/ eff , e^) = 
(/ eff /|/ eff |)A(|/ eff |,ei), i.e. the phase of the input and 
output are directly related, although not equal. 



A. Summary of clean JCHM 

For comparison, we briefly summarize here the behav- 
ior in the absence of disorder. In the absence of hopping, 
the problem is identical to that studied by Bishop et al 
(34j : The coupled qubit-cavity system has an anharmonic 
polariton spectrum, and so at low pumping, one can con- 
sider the response to pumping an effective two-level sys- 
tem. If one considers the coherent field amplitude | (a) | as 
a function of pump frequency u p then at weak pumping 
there is a standard Lorentzian response, while at higher 
power, power broadening [35| leads to a reduction of the 
coherent field amplitude near resonance, i.e. there is an 
anti-resonance feature. Turning on hopping, the location 
of the anti-resonance shifts away from the low-power res- 
onance. Eventually it shifts so far that the coherent field 
amplitude vs pump frequency develops a jump and an as- 
sociated bistability. Such bistability is analogous to that 
known in the Dicke model when driving above resonance, 
where nonlinearity can blueshift the polariton frequency 
into resonance. 

Let us note at this stage that although the existence of 
bistability is due to the mean-field decoupling, its pres- 
ence is indicative of physically meanin gful bimodal distri- 
butions in the true density matrix [36|, [37j . The equation 
of motion for the full-system density matrix is linear, and 
so either has a unique steady state, or a degenerate sub- 
space of steady states. The mean-field decoupling instead 
produces a nonlinear equation for the single-site den- 
sity matrix, which may have multiple distinct solutions 
— these distinct solutions can thus describe bistability. 
Where mean-field theory would predict bistability, the 
full density matrix would generally have a configuration 
with a significant weight near both of these mean-field so- 
lutions, but with a fixed ratio between their weights and 
a tail of finite probability states that connect these. Both 
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the ratio of weights and the existence of the intermedi- 
ate states cannot be found by mean field theories, and 
require consideration of fluctuations, and specifically in- 
stanton and soliton corrections that would describe tun- 
neling between different mean-field configurations [38| . It 
is however worth noting that all these statements relate 
to the ensemble averaged steady state of the system. If a 
system is prepared near to one of the two bistable states, 
the subsequent dynamics will initially remain near that 
configuration until a tunneling event causes a transition 
to the other state. Such tunneling (quantum, thermal or 
induced by external noise) can cause transitions in both 
directions, and eventually produces a fixed ratio between 
the two parts of the bimodal distribution. 

Since the spacing of energy levels of the JCHM is an- 
harmonic, in the limit of relatively weak hopping, the 
problem can be understood quantitatively by restricting 
the on-site problem to a reduced Hilbert space of 0, 1 ex- 
citations. As discussed in [28|, |34[, this is valid as long 
as other excitations are sufficiently far from resonance, 
U e ff ^> f where U e s is an effective anharmonicity (which 
vanishes for large hopping). This reduces the problem 
to: 



eff 



-U 



v 



M 



(ij) 



Z 



(10) 



where rf are Pauli matrices in the reduced Hilbert 
space and the effective parameters are r] = (J + e — 

J = J sin 2 0, and / = — / sin # 
e). Losses are described by 



v /(J-e) 2 + 4^)/2-o; p 
with tan(26>) = 2g/(J - 

Z) i («/2)£[r i "] with k = ft sin 2 + 7 cos 2 (9. Since J <C g 
is assumed one may further approximate 77 c± — g + ( J + 
e)/2 — ~ 7r/4. The steady state of this problem can 
be reduced to coupled equations for the coherent field am- 
plitude ip = (r~), an effective detuning A = 77+ J(2n— 1), 
and the excited state population n = (1 + t 2 )/2, 



/(A - in/2) 
A 2 + (k/2) 2 + 2/ 2 



P 



A 2 + (k/2) 2 + 2/ 2 



( n ) 

One may thus see that for / <C k one has resonance at 
r] = J = J/2 giving uj p c± —g. In contrast, for larger / 
one has n 1/2 and the center of the anti-resonance is at 
77 = 0, i.e. Up = —g + J/2. Such behavior is already clear 
in the clean limit shown in Fig. [TJ even with / = k = 7. 

For large enough J, there are multiple solutions of the 
above equations. Equivalent ly this means 77 is a non- 
monotonic function of A, and so by writing: 



77 = A + J - J 



2f 



A 2 + (k/2) 2 + 2f 2 



(12) 



one can find the critical value of J for bistability by seek- 
ing the solution of dr]/dA = = d 2 r]/dA 2 . This yields: 



(13) 




For J > J c , there is a range of 77 (i.e. pump frequencies) 
for which A (77) and thus ^(rj) are multi- valued and so 
describe bistability. 



B. Effects of disorder 

We consider first the effects of disorder when J < J c , so 
there is no bistability, but a strong distortion compared 
to J = 0. The probability distribution of the amplitude 
of the coherent field strength, P(\ip\) in this case is shown 
in Fig. [TJ and cross sections, showing the full probability 
distribution P(ip) over the complex plane are given in 
Fig. [21 We consider here parameters as discussed in [28| 
for ease of comparison. For the inclusion of disorder, 
one may note that a typical scale of disorder in recent 
experiments [39| is <j € ~ 1 MHz, corresponding to 0.002 < 
Ve/g ^ 0.005. We show results for a e /g ~ 0.002; larger 
disorders show very similar behavior. 
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FIG. 1. (Color online) Probability distribution of as a 
function of pump frequency. The dashed (cyan) line indicates 
the value of in the clean limit, and the colormap shows 
the probability distribution of for Gaussian disorder of 
variance a € / 'g = 0.002g. The dash-dotted (gray) line indicates 
the approximate solution to the clean case given above in 
Eq. ©. Blue arrows mark the values of pump frequency at 
which the full probability distribution of complex ip is shown 
in Fig. [21 Other parameters are / = k = 7 = 0.005g, J / g — 
0.020 and a geometry with z — 2 is assumed. 

For most pump frequencies, disorder has a relatively 
weak effect, but near the anti-resonance feature it causes 
a much larger effect. This can easily be understood from 
the discussion of the clean case above: in this regime 
uj p ~ — g + J/2, and the effective detuning ~ rji c± 
—g J r(J J r€i)/2 — ujp ~ Ci/2. Thus, near the antireso- 
nance, the variance of A is large compared to its mean 
value. Since the variance of disorder is of the same order 
as the linewidth ft, one finds in this regime that the phase 
of the on-site order parameter can vary significantly. This 
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is clearly seen in Fig. Efb). In contrast, away from this 
point, the mean value of A is much larger than its vari- 
ance, and so disorder has only a weak effect on the phase 
and amplitude, hence the clean results are recovered. 
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FIG. 2. (Color online) Probability distribution of the complex 
observable ip at four different pump frequencies for J < J c . 
The blue crosses indicate the value of ip found for the clean 
case, and the colormap shows the probability distribution 
with (T e /g — 0.002. All parameters are as in Fig. [1] 

As one continues to increase the pump frequency above 
the anti-resonance, the field amplitude remains notably 
higher than in the clean case, and (as seen in Fig. [2jc)) 
the phase distribution remains broad. The increased 
amplitude can be clearly understood as an effect of 
the phase distribution: increasing the phase distribution 
means the convolution distribution Q(</>) moves toward 
smaller <j). Since the field seen by a given site is given 
by / eff = f — J4>/z, and since Re((/>) > 0, reducing \(j>\ 
increases the effective driving, and thus increases the am- 
plitude. 

The phase spreading seen here signifies an important 
distinction between the thermal and the non-equilibrium 
disordered problem. In the thermal case, a real distribu- 
tion of ip is stable, but in the non-equilibrium case there 
is always a distribution of phase, and near resonance, this 
becomes particularly notable. In the current case, phase 
symmetry is broken by the external pump. However, 
for incoherent pumping, phase symmetry is not broken. 
The presence of phase spreading then means that follow- 
ing Imry and Ma [40|, no spontaneous phase symmetry 
breaking is possible in d < 4. This is quite different from 
the equilibrium JCHM where a superfluid (superradiant) 
state with phase symmetry breaking is expected in d > 2. 
A similar observation has recently been made for the dis- 
ordered polariton condensate [41( . 

As discussed above, in the clean case, for J > J c bista- 
bility occurs because of the multivalued nature of A (77). 
However, the range of detunings where this occurs is the 
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FIG. 3. (Color online) Probability distribution of as a 
function of pump frequency for J > J c . Lines and parameters 
are as for Fig. [TJ except J j g — 0.04 in this case. 



same range where strong phase spreading was seen, and 
thus disorder strongly affects the behavior in exactly this 
region. Thus, as seen in Fig.[3j the disordered case with a 
typical disorder strength cr e /g = 0.002 does not show any 
bistability, with no weight near the new clean solution 
which appears as uo p is increased. The absence of bista- 
bility is revealed by noting that the same steady state 
is found independent of starting distribution of P($)\ in 
the current case this was tested by comparing a "sweep" 
of slowly increasing or decreasing pump frequency, both 
cases lead to identical results. Since a unique solution 
exists in these cases, this may be taken as an approx- 
imation of the true asymptotic state of the disordered 
system. One may also note in Figs. [4f e-h) that even 
once the low uo p solution vanishes, weight is not concen- 
trated near the high uj p solution until significantly above 
the antiresonance frequency, as there is a strong effect of 
the phase distribution. 



IV. CONCLUSIONS AND OPEN QUESTIONS 

We have presented a non-equilibrium extension of 
stochastic mean-field theory, applicable to problems of 
coupled cavities with rather general forms of driving and 
dissipation. Using this approach we studied the effect of 
disorder on the driven dissipative JCHM. Near the anti- 
resonance, disorder introduces significant phase spread- 
ing, which in turn increases the coherent field ampli- 
tude over a range of pump frequencies above the anti- 
resonance. 

The results presented above for the dissipative driven 
JCHM clearly demonstrate that the combination of open 
quantum systems with disorder can lead to behavior that 
is not seen with only one of these two ingredients in iso- 
lation. Such behavior prompts an important question 
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FIG. 4. (Color online) Probability distribution of the complex 
observable ip at nine different pump frequencies for J > J c . 
All parameters are as in Fig. [3] 



regarding whether dissipative coupled matter-light sys- 
tems could ever be used as "quantum simulators" of dis- 
ordered models. At the same time, it indicates that there 
are open questions as to what the phase diagram of inco- 
herently pumped disordered dissipative systems may be. 
In some cases, such as the Rabi-Hubbard model [26], only 
discrete symmetries exist, and so the effects of disorder 
should not destroy the symmetry-broken phase, and the 
behavior may be equivalent to that of the site-disordered 
transverse field Ising model. However, for cases with 
continuous symmetry it is unclear whether any phase 
boundaries exist, since neither the superfluid nor Mott- 
insulating states survive the effects of dissipation. Such 
questions can be in part addressed by the SMFT ap- 
proach described here. 

Another challenge is to go beyond the Stochastic- 
Mean-Field limit presented here, and produce alternate 



methods to treat open disordered lattice problems. As 
with all mean- field approaches, SMFT neglects quan- 
tum correlations between different sites; an assumption 
only valid in the limit of high coordination. In addition, 
SMFT makes a second assumption, that there is no cor- 
relation between the on-site energy and the field distri- 
bution seen. Such an assumption implies self-averaging, 
while it is known that self- averaging breaks down in the 
equilibrium Bose glass [23]. An alternative approach that 
may circumvent this is to consider extensions of the cav- 
ity method e.g. [Hj. For the purpose of understanding 
the behavior of currently achievable experiments Q, fi- 
nite size simulations of the mean-field [28| or beyond- 
mean- field [lil, [30| dynamics may be more appropriate. 
However, a full understanding of the behavior of such 
dissipative models may well depend on rare events, not 
captured in finite size simulations, so methods such as 
that presented here may play an important role. 

In conclusion, the combination of dissipation and dis- 
order can lead to types of behavior in coupled cav- 
ity arrays that cannot be seen in either the clean non- 
equilibrium system, or disordered equilibrium case. This 
suggests that such cavity arrays may not be appropriate 
as quantum simulators to understand equilibrium disor- 
dered problems. Stochastic mean-field theory can pro- 
vide a simple route to address some classes of system, 
but leads to questions that require more sophisticated 
approaches to non-equilibrium disordered problems. 
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